home *** CD-ROM | disk | FTP | other *** search
- /*
- * Code to generate bezier control points to interpolate between points.
- * Copyright Shamim Mohamed 1991
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
- *
- * shamim@cs.arizona.edu
- * Shamim Mohamed
- * Dept of Computer Science
- * University of Arizona
- * Tucson, AZ 85721
- *
- * --------------------------------------------------------------------------
- *
- * Inputs: an array of points, the number of points, two fudge-factors,
- * and two procedures that are called for moveto and curveto. I use d = 3.0
- * and e = 0.2. d controls the slope of the curve depending on the neighbours
- * (it's not prob. worth futzing with it) and e controls the shape of the ends
- * of the curve. Experimenting with e may be fruitful.
- * This is NOT an approximation algorithm. Given n points, it ALWAYS generates
- * n-1 splines that pass through all the given points. For some sort of
- * optimum approximation, see the paper in "Graphics Gems", Klossner ed. (It's
- * the last source-code in the Appendix.)
- *
- * The callback procedures can be something like:
- * void moveto(x, y)
- * double x, y;
- * {
- * printf("%g %g moveto\n", x, y);
- * }
- * void curveto(x1, y1, x2, y2, x3, y3)
- * double x, y;
- * {
- * printf("%g %g %g %g %g curveto\n", x1, y1, x2, y2, x3, y3);
- * }
- *
- */
-
- #include "plotfoil.h"
- #include "externs.h"
-
- void interpolate(p, npoints, d, e, moveto, curveto)
- point_t *p;
- int npoints;
- double d, e;
- void (*moveto)(), (*curveto)();
- {
- int i;
- point_t slope_prev, slope_next;
- point_t cp[MAXPOINTS], cpp[MAXPOINTS];
- double l[MAXPOINTS];
-
- l[1] = sqrt(sqr(p[0].x-p[1].x)+sqr(p[0].y-p[1].y));
- slope_prev.x = p[1].x - p[0].x;
- slope_prev.y = p[1].y - p[0].y;
- for(i = 1; i < npoints-1; i++) {
- point_t slope;
- double m;
-
- l[i+1] = sqrt(sqr(p[i].x-p[i+1].x)+sqr(p[i].y-p[i+1].y));
-
- slope_next.x = p[i+1].x - p[i].x;
- slope_next.y = p[i+1].y - p[i].y;
-
- /*
- * Slope of the curve at this point is the weighted ave. of
- * straight-line slopes to neighbours (weighted by length of chord)
- */
-
- slope.x = slope_prev.x * l[i+1] + slope_next.x * l[i];
- slope.y = slope_prev.y * l[i+1] + slope_next.y * l[i];
-
- slope_prev = slope_next;
-
- /*
- * Alternative method:
- * Parameters:
- * a 0 Higher number makes curve tighter. -1<a<1.
- * b .1 Higher number makes sharp curves tighter
- * and broad curves less tight. -1<b<1
- *
- * Given points A B and C, bezier control points for B are aligned along
- * a line that is perpendicular to the line the bisects the angle ABC.
- * The distance of the control points from B is determined by the
- * following formula that was found by trial and error is has no
- * particular theory behind it:
- * Distance of control point from B towards A =
- * (1/2)*(distance between A and B)*
- * (cos(ABC)^(a*b + a + b))(b + 2(1-b)(1-a)/(3-a))
- * (Timothy van Zandt <tvz@princeton.edu>)
- * My opinion: one has to take an inverse cosine and raise it to a
- * fractional power, and the results are very close to my simpler
- * algorithm, so I don't think it's worth it. If someone decides to
- * implement this, please let me know!
- */
-
- /* Now we have slope: m is for the normalisation. */
- m = sqrt(sqr(slope.x)+sqr(slope.y));
- cp[i].x = p[i].x - l[i]*slope.x/(d*m);
- cp[i].y = p[i].y - l[i]*slope.y/(d*m);
- cpp[i].x = p[i].x + l[i+1]*slope.x/(d*m);
- cpp[i].y = p[i].y + l[i+1]*slope.y/(d*m);
- }
- cpp[0].x = (e*cp[1].x + p[0].x)/(1+e);
- cpp[0].y = (e*cp[1].y + p[0].y)/(1+e);
- cp[npoints-1].x = (e*cpp[npoints-2].x + p[npoints-1].x)/(1+e);
- cp[npoints-1].y = (e*cpp[npoints-2].y + p[npoints-1].y)/(1+e);
-
- (*moveto)(p[0].x, p[0].y);
- for(i = 1; i < npoints; i++)
- (*curveto)(cpp[i-1].x, cpp[i-1].y, cp[i].x, cp[i].y, p[i].x, p[i].y);
- }
-
-